Synchronization of Chaos in Fully-Developed Turbulence 
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We investigate chaos synchronization of small-scale motions in the three-dimensional turbulent 
energy cascade, via pseudo-spectral simulations of the incompressible Navier-Stokes equations. The 
modes of the turbulent velocity field below about 20 Kolmogorov dissipation lengths are found to 
be slaved to the chaotic dynamics of larger-scale modes. The dynamics of all dissipation-range 
modes can be recovered to full numerical precision by solving small-scale dynamical equations with 
the given large-scale solution as an input, regardless of initial condition. The synchronization rate 
exponent scales with the Kolmogorov dissipation time-scale, with possible weak corrections due to 
intermittency. Our results suggest that all sub- Kolmogorov length modes should be fully recoverable 
from numerical simulations with standard, Kolmogorov- length grid resolutions. 



Chaos synchronization (CS) [1] is an intriguing phe- 
nomenon which has been defined as "a process wherein 
two (or many) chaotic systems . . . adjust a given prop- 
erty of their motion to a common behavior due to a cou- 
pling or to a forcing" [2]. The simplest example is a 
chaotic dynamics x = f(x) whose phase vector x is pro- 
jected onto two orthogonal components xi = Pix and 
x i = Qix = x — xi satisfying two coupled equations 



^x 1 = P 1 f(x 1 +x' 1 ), 
^x' 1 = Q 1 f(x 1 +x' 1 ). 



(1) 



Chaos implies sensitive dependence to initial data, with 
nearby trajectories diverging exponentially. However, 
consider another dynamical system in the Qi~ s P ace given 
by a copy of the second equation: 



-w = Q 1 f(x 1 (t)+w) 



(2) 



with Xi(i) substituted from the solution of (1). Chaos 
synchronization occurs if the trajectories w(i) and x i(t) 
converge, limt-^oo ||w(t) — x^(t)|| = 0, for an arbitrary 
choice of initial condition wo in (2). Such a phenomenon 
requires that the leading Lyapunov exponent for the sub- 
dynamics (2) be negative. It is often the case that syn- 
chronization occurs, at least approximately, even when 
imperfect data xi(£) is employed in (2), e.g. the exact 
xi(i) contaminated with substantial errors. This effect 
was proposed in [3-5] as a basis for encrypted commu- 
nications. CS has also been observed in neural networks 
[6, 7], with "hyper-synchronous" dynamics in the human 
brain associated to epileptic seizures [8]. CS has been re- 
ported in spatio-temporal chaos [2, 9-11], investigated for 
fluid models used in meteorology [12-14], and proposed 
as a mechanism for turbulence control [15-17]. 

No previous numerical study of CS has been made, 
to our knowledge, for fully-developed three-dimensional 
Navier-Stokes (NS) turbulence with a Kolmogorov iner- 
tial range. Our goal in this Letter is to explore CS for 
3D NS with Pi taken to be the projection onto the fi- 



nite number of velocity modes with wavenumber magni- 
tudes less than a fraction / of tt/tjk, where t\k is the 
Kolmogorov disspation scale, and with Qi the orthog- 
onal projection onto the modes with higher wavenum- 
bcrs. Our principal motivation is experimental results 
[18], theoretical work [19, 20], and numerical simulations 
[21, 22] implying that spatial intermittency can lead to 
length scales far smaller than the Kolmogorov scale tjk- 
It has been argued on the basis of such tiny unresolved 
length scales "that the DNS [direct numerical simula- 
tion] based on the mesh equal to the Kolmogorov scale 
becomes quite inaccurate" [20]. If true, this would call 
into question the vast majority of current DNS studies 
of turbulent flow. A contrary argument is based on the 
idea that the sub-Kolmogorov scales should be "slaved" 
to the inertial-range modes and, thus, implicit and recov- 
erable from DNS with grid resolution tjk- A mathemati- 
cal formalization of this idea closely related to CS is the 
notion of an inertial manifold (IM) [23], which consists 
of an invariant, attractive manifold given by the graph 
of a mapping x^ = <&(xi) which recovers x^ for given 
Xi. Existence of an IM with the property of "asymptotic 
completeness" [24] is one possible mechanism for CS (e.g. 
see [25]). There are currently no proofs of existence of 
an IM for 3D NS dynamics, although "approximate IM" 
have been obtained for 2D NS [26, 27]. These have been 
proposed for use as nonlinear Galerkin approximations 
to the dynamics of "large" super-Kolmogorov scales in 
NS turbulence, whereas our goal is the opposite one to 
recover the sub-Kolmogorov scales and address the out- 
standing issue of the smallest length-scale in a turbulent 
flow [21, 22]. Even when existence of an approximate IM 
can be established for NS, there are no sharp estimates of 
the smallest fraction / of the Kolmogorov wavenumber 
sufficient for slaving. Thus our numerical investigation 
in this Letter is an important complement to existing 
mathematical results. 



The incompressible NS equations with a solenoidal 
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body force f have the form: 



d t u + P(u • Vu - vAu) = f , 



(3) 



where u is the fluid velocity, v is the kinematic viscosity, 
and P is the Leray projection to enforce the incompress- 
ibility condition V • u = 0. We employ the pseudospectral 
DNS method, which solves a Galerkin approximation to 
this system 



d t u 2 + -P2(u 2 • Vu 2 - i/Au 2 ) = f. 



(4) 



with Pi the Leray projection in a space spanned by a 
finite set P>i of Fourier modes. For our purpose, this will 
represent the "fine-grained" NS dynamics. In order to 
study the possible presence of CS, we consider a further 
subset P>\ C i?2 and corresponding projections P\ and 
Qi = P% — P\. The subdynamics for the evolution of 
w = U2 — Ui is given by the equation: 

d t w + Qi [(ui + w) • V(ui + w) - i/Aw] = 0, (5) 

where we have assumed that f forces only the large scales, 
i.e. Qif = 0. In our experiments we shall solve the fine- 
grained equation (4) for 112(2) and then solve the sub- 
dynamics (5) with Ui(t) = Pi 112(2). We shall investigate 
whether w(2) converges to Qiu 2 (2) for increasing 2, inde- 
pendent of the initial data w . Specifically, we will study 
the evolution of the normalized error defined as 



e(2) 



|w(*)-Qni 2 (t)|j 2 
l|w(2)|| 2 



(6) 



where ||||2 is the usual L 2 norm. If synchronization oc- 
curs, this error should tend to zero exponentially fast, 
independent of the initial value Wo. 

The concrete system considered in this work is Kol- 
mogorov flow with f = [Asm(kfy), 0, 0) for A = 1 and 
kf = 1, in an elongated box [0, L x ] x [0, L y ] x [0, L z ] with 
L x = 3L y = 3L Z = 6tt. The numerical simulation uses a 
space grid of N = 3n x n x n points with isotropic mesh- 
spacing corresponding to maximum wavenumber fc^/ = 
n/2. This particular configuration leads to a nontrivial 
turbulent flow, that is anisotropic and inhomogeneous in 
the large scales [28]. In the context of synchronization of 
chaos, it is relevant that strong bursts can be observed 
in Kolmogorov flow. In [28] very long integration times 
were used precisely because the time averages presented 
converge very slowly. Thus by studying the system at dif- 
ferent times, significantly different regimes can be sam- 
pled. In terms of the kinetic energy E and energy dissipa- 
tion rate e, the Kolmogorov units and the Reynolds num- 



ber are r\ K = (^f\ 



IE 



3 (wr)i/2 

Five series of simulations are performed, with resolutions 
ranging from 144 x 48 x 48 to 768 x 256 x 256 grid points, 
and R\ going from 40 up to 250, keeping the minimum 
^mVk around 1.5. The energy spectra plotted in Fig. 1 
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FIG. 1. Energy spectra of the Kolmogorov flow simulations 
for n — 48, 64, 96, 128, 256. The spectra are taken from in- 
stantaneous snapshots with no time-averaging. The moderate 
quality of the collapse in dissipation-scale units is likely due 
mostly to strong unsteadiness in the flow. 



show a short Kolmogorov inertial range with approximate 
—5/3 power- law scaling. 

For our CS study a very long simulation of Kolmogorov 
flow is performed for each resolution, saving a few time 
series of the velocity fields from the quasi-stationary 
regime, each interval separated by relatively long times. 
Four time intervals of 112 are chosen for each resolution 
(three for the 768 x 256 x 256 case). Next Ui is obtained 
by the projection Pi of u 2 onto modes with wavenum- 
bers smaller than a cutoff value in each direction (i.e. 
\k x \, \k y \, \k z \ < k c ). Finally, w is evolved in time using 
(5). For each interval, two initial conditions Wo were cho- 
sen, so that each series consists of eight individual runs. 
In the experiments presented, initial data wo with "natu- 
ral" spectral scaling properties were created by applying 
random phase shifts to all Fourier modes of 112 ui. 
Several alternative initialization methods for wo were 
tested and yielded consistent results, not shown here. As 
observed in Fig. 2, for the indicated values of k c , e(2) 
does indeed decrease exponentially fast, until it reaches 
a smallest possible value dictated by our single precision 
arithmetic. Thus w synchronizes to Qiu 2 . 

Fig. 2 also shows that the exponential decay rate a 
becomes greater at larger k c , a natural result since w 
then lives on smaller and hence faster scales. We have 
studied this effect quantitatively The linear part of the 
trends in Fig. 2 can be computed from the data by least- 
square error fitting at + b to the measured log 10 e(2) in 
the region where the error is larger than the roundoff 
error floor, i.e. for e between 1 and 10~ 5 . The behav- 
ior of the measured a as function of k c depended on the 
various parameters of the simulations. To attempt to 
collapse the results, various non-dimensionalizations for 
a and k c were tested. It was found that good collapse 
is observed when using Kolmogorov (viscous) scales for 
both the cutoff wavenumber as well as the synchroniza- 
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FIG. 2. Time evolution of normalized error e(t) for the simu- 
lation on a grid of 308 x 96 x 96 and R\ ~ 108, using several 
cutoff-wavenumbers. The slope in these graphs yields the ex- 
ponential decay rate a, or rate of synchonization. 



tion exponent, i.e. to plot arx versus k c r\K- See Fig. 3. 
To document the scatter due to possible lack of statisti- 
cal convergence, the duration of "exponential decay" was 
split in half for each individual run, the corresponding 
pair (/c c ?7, o,tk) was computed for each of the resulting 
e(t) histories, and then the average over all the simu- 
lations with the same resolution (or Reynolds number) 
was computed. These are the results that are presented 
as symbols in Fig. 3. Error bars are for maximum and 
minimum values. The results collapse reasonably well, 
although the lines seem to shift a little to the right with 
increasing resolutions. This hints at a slight Reynolds 
number dependence, which is expected due to intermit- 
tency [19, 20]. The results of Fig. 3 are parameterized 
well by a linear fit ar^ ~ — P{k c r}K~ 0.15), implying that 
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FIG. 3. Symbols: average of measured synchronization ex- 
ponents (obtained by fitting the exponential range in results 
such as in Fig. 2) as function of cutoff-wavenumber for five 
different simulation sizes and Reynolds numbers, plotted in 
Kolmogorov units. Error bars are for maximum and mini- 
mum values of different runs and ranges used in the fit. 




FIG. 4. Kinetic energies, taken for a fixed z from a 768 x 
256 x 256 simulation (x varies on the vertical and y on the 
horizontal). Left: coarse grained field obtained with k c r/K ~ 
1/4. Right: refined version of coarse grained field. The values 
are normalized with the volume averaged kinetic energy of 
the original field. Note that these snapshots were taken after 
synchronization had taken place, so ui + w is equal within 
numerical precision to the original field U2 . The latter cannot 
be distinguished by eye from the reconstructed field. 



synchronization of small scales to large scales occurs only 
if the cutoff wavenumber is such that k c r\K > 0.15. Using 
the correspondence r c — this denotes scales smaller 
than r c < 20j]k, i.e. in the transition zone between the 
inertial and viscous ranges. 

The key point to be taken from this study is that it 
is possible to reconstruct perfectly the small scales of a 
turbulent flow from coarse-grained data. If the velocity 
of a turbulent fluid is sampled on a spatial grid even as 
coarse as 10-15 times the Kolmogorov scale, these (time- 
dependent) data can be refined to their original resolu- 
tion, in the sense that the subdynamics of small scales 
after a suitable time will synchronize with the large-scale 
dynamics. Figure 4 shows what this refinement implies: 
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fine details of small-scale structures that are smeared out 
in the coarse-grained field reappear, as if by magic, when 
refined by computing the subdynamics. Of course, syn- 
chronization takes time. For example, assuming that 
kcVK ~ 1/4 as in Fig. 4, and assuming that a preci- 
sion of e = 1CT 3 is desired, an interval of about 50tk is 
needed. This translates into about 65/ R\ in units of the 
integral time, significantly less then an integral time for 
moderately large values of R\. Doubling of this interval 
would lead to an error of order 10 -6 , at the lower limit 
for single precision computations. 

Our results offer some support to the current prac- 
tice of DNS with grid spacing of order t\k, since they 
suggest that there may be an exact solution of 3D NS 
which, when coarse-grained to the grid scale, agrees 
with the finite-resolution simulation. Tiny scales much 
smaller than t]k may be present, but completely slaved 
to the super-Kolmogorov scales. To more fully address 
these issues, numerical experiments on CS must be per- 
formed with approximate data for Ui which come not 
from a projection of a fine-grained solution U2 but in- 
stead from a pscudospcctral DNS with cutoff wavenum- 
ber k c . Outstanding issues arc whether CS will occur for 
such approximate Ui and whether the reconstructed field 
U2 = Ui + w is then a solution of the fine-grained equa- 
tions. These questions are currently under active investi- 
gation. The size of the smallest length-scale in turbulence 
is of interest not only for physical theory but also for fun- 
damental mathematical theory of 3D incompressible NS. 
The Clay Millenium Prize problem on that equation con- 
cerns whether its solutions at sufficiently high Reynolds 
numbers may develop actual singularities, with veloci- 
ties exploding to infinity at the singularity and smallest 
length scale going to zero [29]. In nature, physical effects 
beyond incompressible NS would cut off the singularity at 
some tiny length-scale, but the observable manifestations 
should be striking. There is presently no empirical evi- 
dence whatsoever for such "Leray singularities", but this 
may be due to limited resolution or statistics of current 
numerical and experimental studies. If such singularities 
occur anywhere at all, high Reynolds turbulent solutions 
are perhaps the most likely venue. Better understanding 
of the interactions between inertial range and far dissi- 
pation range modes in turbulent NS flows should help to 
illuminate this problem. 
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